From signal-based to comprehensive magnetic resonance imaging

We present and evaluate a new insight into magnetic resonance imaging (MRI). It is based on the algebraic description of the magnetization during the transient response—including intrinsic magnetic resonance parameters such as longitudinal and transverse relaxation times (T1, T2) and proton density (PD) and experimental conditions such as radiofrequency field (B1) and constant/homogeneous magnetic field (B0) from associated scanners. We exploit the correspondence among three different elements: the signal evolution as a result of a repetitive sequence of blocks of radiofrequency excitation pulses and encoding gradients, the continuous Bloch equations and the mathematical description of a sequence as a linear system. This approach simultaneously provides, in a single measurement, all quantitative parameters of interest as well as associated system imperfections. Finally, we demonstrate the in-vivo applicability of the new concept on a clinical MRI scanner.

the signals are matched to an a priori dictionary to find the relaxation properties. This matching procedure avoids the fitting to an explicit signal model based on repeatedly solving the Bloch equations. It thereby essentially extends the successful compressed sensing principle 16 to the temporal direction.
Some improved alternatives have recently been published, such as Quantitative Transient Imaging (QTI) 10 that uses a fixed time of repetition and a linearly increasing variable flip angle (vFA) and Magnetic Resonance Field Mapping 17 which additionally provides B 1 and B 0 maps through including multiple B 0 and B 1 . values in the dictionaries. Sbrizzi et al. recently proposed a brute force approach named MR-STAT providing all the relevant maps (T 1 , T 2 , PD, B 1 and B 0 ) that also avoids dictionary matching 18 .
In this work, we propose a new MR method based on a mathematical closed-form description of the magnetization evolution of a single species along a sequence of repetitive blocks which contain RF pulses and readout gradients. This provides four important advantages. First, it enables a comprehensive simultaneous estimation of all intrinsic parameters and avoids confounding by experimental imperfections such as B 0 inhomogeneities and B 1 inaccuracies. Second, the analytical description of the sequence of blocks gives a new insight that facilitates the selection of pulse sequence design (TR, FA and phase) to increase sensitivity for all relevant parameters. Third, no dictionary is required to perform the estimation. Fourth, we also show that, the proposed design, based on the insight given by the algebraic mathematical model, avoids the banding artefact-a long-standing issue in the standard balanced Steady State Free Precession (bSSFP) MRI technique.
We propose that the use of the comprehensive MR method described in this work could contribute to standardized MR through the adoption of multiparametric methods in clinical protocols and to increased reproducibility in follow-up and multi-site or multi-vendor studies.

Theoretical framework-from linear algebra to dissipative coupled harmonic oscillators
We provide an algebraic description of the signal evolution during the entire MR sequence. Our goal is to establish a theoretical framework for quantitative sequences that utilizes a single transient response and can yield quantitative maps for all relevant parameters, intrinsic ( T 1 , T 2 , PD) and experimental ( B 0 , B 1 ) at once. First, we will consider an imaging voxel as homogeneous and represented by a single magnetization vector, i.e. we use a single species model. In general, the assumption of the intra-voxel homogeneity is not valid. However, as we demonstrate, it can be maintained by a careful choice of the acquisition scheme.
Discrete algebraic description of an MR imaging sequence. MRI sequences consist of a train of acquisition blocks (Fig. 1). The acquisition blocks alter the magnetization vector. Each block consists of radiofrequency excitations (Fig. 1a, in red), magnetic field gradients (Fig. 1a, in blue) and wait times (Fig. 1a, µperiods for free precession, absent of RF excitations or magnetic field gradients) each of which can be described as linear operator that acts on the magnetization vector 5 . Hence, the combined effect is also a linear operator called propagator.
We focus on a special case when the propagator, A , is identical along the train of acquisition blocks and has duration T R . As only the net effect of the propagator is equal, blocks differing in spatial coding k-space trajectories are possible. In this case the magnetization vector m evolves according to the recursive equation: where m ss is the steady state magnetization and m 0 the initial state. We can focus on the homogeneous recursive Eqs. (5,9): with the solution: where µ n = m n − m ss .
In comparison, for continuous modeling, the kinetics are described by the differential equation: µ n = A n · µ 0 Figure 1. (a) Illustration of the "propagator": a single block of events. μ is the input-output magnetization, in red the RF excitations and in blue the gradients. (b) A typical train of blocks in a MR pulse sequence. We will use identical blocks or operators A = A 1 = A 2 = … = A n .  Figure 2 depicts both the continuous and the discrete solutions.
In the case of an imaging sequence consisting of k consecutive excitations and free precession periods, A takes the following form: with R s (ϑ) specifying a rotation of angle ϑ around axis s ( x , y or z in Eq. 5), where the α, β and γ are a function of the RF flip angle, phase and local off-resonance frequency, and represents the relaxation in the period τ j , j τ j = T R .
A real-valued expression for the signal in an MR sequence. In order to exploit the information content of the signal evolution, it is convenient to have a closed-form expression for µ n in Eq. (3). This requires an expression for A n , where A is a general, real-valued 3 × 3 matrix.
The usual method to derive an expression for A n is to diagonalize A by its eigenvectors provides a simple expression with the eigenvalues 1 , 2 , 3 : (4) d dt y(t) = B · y(t), withthesolution : y(t) = A t/T R · y(0)) (a-c) The repeated block consists of an excitation with flip angle α=30º, β=14º (accumulated phase as a result of off-resonance during T R ,T R =10 ms), relaxation times are T 1 = 878 ms, T 2 = 47.5 ms. n 1 and n 2 span the plane of oscillation. The magnetization vector always points to a point of this plane. Only the orientation of the plane is fixed throughout the evolution; it shifts parallel to n 3 .   3 ] is the matrix formed by the eigenvectors of A 5,8 . A sufficient and necessary condition of the diagonalization is det(V ) = 0 , i.e. the eigenvectors of A are linearly independent. This condition is met if the eigenvalues are distinct, however this is not a necessary condition (multiplicity of the eigenvalues of A alone does not preclude diagonalization).
We follow a different method for derivation of the expression of A n . This method does not rely on the diagonalizability of A and it is also idependent of the multiplicity of the eigenvalues. The discrete Putzer's algorithm 20,21 relies on a consequence of the Cayley-Hamilton theorem: A n can be expressed as an ( d-1) order polynomial of A ∈ R d×d , in our case with the identity matrix I , A and A 2 . Applying the Putzer's algorithm the resulting expression for 3 × 3 matrices with distinct eigenvalues is: The expression is also valid in its limit for multiple eigenvalues, may any of the eigenvalues as roots of the characteristic polynomial be identical: A is a product of rotation and dilation operations, represented by a non-singular, real valued matrix. The eigenvalues of A can always be expressed as 1 = ρe iϕ , 2 = ρe −iϕ , and 3 = η , where ρ and η are real valued, and ϕ is either real or purely imaginary. A has one real-valued eigenvalue and its real eigenvector ( 3 and v 3 ), but the other two eigenvalue and eigenvector pairs are usually complex conjugates. It is also relevant to realize that ρ 2 · η is the determinant of the propagator which depends on the intrinsic T1 and T2 and it is essentially the total dissipation during the entire propagator block. This amount can be voxel-wise obtained and it can be reconstructed as a map. This map is characteristic of the specific harmonic oscillator and we will refer to it as the HO (harmonic oscillator) map, or image, for the rest of the document.
Substituting the eigenvalues in Eq. (8), with a somewhat lengthy, but otherwise straightforward derivation we get: The real-valued matrices M j are: The expression for the special cases when two or all three eigenvalues are identical can be derived in the limit, e.g. ϕ → 0 , ρe iϕ → η.
So, the explicit discrete expression of the magnetization evolution as in Eq. (3) is : where the real-valued normal mode 22 vectors n 1 = M 1 · µ 0 , n 2 = M 2 · µ 0 and n 3 = M 3 · µ 0 depend on the initial state µ 0 . At this point it is important to make two important remarks. First, in order to clarify the terminology: these normal modes are not the eigenvectors (sometimes also called normal modes, especially if A is a symmetric matrix). They are always real valued and represent special directions in real 3D space, along which motion is defined in Eq. (11).
Second, it is important to point out that the matrix expressions in Eq. (10). include a projection matrix: . Based on this one, we can predict that the normal modes can also have zero length. This happens to be the case when µ 0 is in the null-space of a matrix M 1 , M 2 or M 3 . The normal modes n 1 , n 2 and n 3 form a complete non-orthogonal basis in 3D if µ 0 is not in the nullspace of any of the matrices in Eq. (10).
It is clear from Eq. (11) that µ follows decaying oscillations in two directions n 1 and n 2 with locked phase and frequency. In the third n 3 direction it will follow an exponential decay towards zero.
The detectable signal is a projection of m n on the plane of detection ( xy-plane). It is a straightforward calculation to show from Eq. (11), by the usual quadrature detection with a complex expression: where a , b , δ and ξ are real parameters determined by the initial value m 0 and the propagator. δ is the phase difference between transmit and receive. This equation resembles Torreys' solution of the Bloch Eq. (3) but stretched to the repetitive pulse sequence.

Experimental design
There are three conceptual aspects of our experimental design that render our technique viable in terms of the applicability of a single species model and the sensitivity to intrinsic parameters: (a) Information content of the signal, (b) Single species description, (c) Eigenvalues of the propagator.
(a) Information content of signal-resonance condition. The description of the motion of µ allows the physical analogy of three independent, linear and damped oscillators, linked to each other only by the initial condition 9,23 . Although n 1 , n 2 and n 3 span the space that holds the trajectory of the magnetization, they do not necessarily form a complete basis in 3D. In Fig. 2, the trajectories of µ are depicted in the rotating frame of resonance, where the propagator consists of a single excitation. The case of off-resonance ( β = 0 ) and on-resonance ( β = 0 ) are shown, where β is the phase accumulated during a TR. In practice, B 0 can not be fully controlled in an imaging volume, and β can take any value in [0, 2π] . In the on-resonance case, the last term in Eq. (11) vanishes because n 3 = 0 , therefore the signal evolution carries no information about the real eigenvalue η . The information loss on-resonance can not be avoided with a block containing only one excitation. Hence, a special composite propagator is required.
(b) Single species model. Equation (11) describes the temporal evolution of the magnetization in a voxel under the assumption that it can be characterized as a single species. There are arguably two reasons why this assumption is not valid: (1) Related to limitations of the imaging technology: Multiple distinct tissue properties in one voxel due to the finite size of voxels or due to experimental imperfections in magnetic and RF fields ( B 0 and B 1 inhomogeneities). (2) Related to the limitations of the Bloch equations, which are an approximation of the complex microscopic behavior of spin relaxation 1 , excluding multiple-pools for magnetization transfer, diffusion etc. In our experimental design we limit the B 1 imperfections with an appropriate slice profile, and the imperfections in B 0 homogeneity on a more conceptual level in the following section.
(c) Eigenvalues of the propagator. An important aspect of adopting an analytically described signal evolution is to estimate intrinsic and experimental parameters. As described in Eqs. (11) and (12), the theoretical description relies on the eigenvalues of the propagator. The propagator can also be viewed as a mapping between the parameter space (T 2 , β) (Fig. 3a) and the space of one of the complex eigenvalues (ρ, ϕ) . As a minimal requirement, this mapping should preserve topology, i.e. it should maintain proximity between points and should be single valued. Otherwise, it would result in indistinguible signal behaviours for distinct (T 2 , β) species and, consequently, the MR experiment would loose sensitivity to some particular combinations of T 2 and β.
Design of a pulse sequence-composite propagator. The propagator that facilitates parametric mapping should fulfill the following requirements: www.nature.com/scientificreports/ • The effect of variation of B 0 on the eigenvalues is limited in order to allow the single species description • The 3D trajectory of µ is maintained at any value of B 0 B 0 in order to avoid loss of information (none of the normal mode vectors vanish) • The mapping between parameter space and eigenvalue space preserves topology in order to distinguish between species by their signal evolution Figure 3 illustrates the problem of non-preserved topology and also demonstrates proper choices of the scheme, which we use for a set of propagators fulfilling the properties for a better estimation. Also, the inspection of the eigenvalue space provides immediate insight into the signal evolution: the radial distance from the origin determines the decay rate, and the angular position determines the oscillation frequency in the signal along the train of acquisition blocks. Figure 3a shows the grid of (T 2 , β) pairs that are evaluated with fixed T 1 and excitation flip angle. Figure 3b-e show the resulting eigenvalue space in the complex plane for four different schemes. Figure 3b shows the " α x " scheme, in which a block consists of a single excitation with 30° flip angle and readout. Figure 3c shows the " α x − α y " scheme that consists of two equally spaced excitations along x and y axis with α = 30 . Figure 3d shows the " α x − α x+δ " scheme consisting of two excitations along the x axis and x + 75 • axis with flip angle α = 150 . Figure 3e shows the " α x − γ y − α y − γ x " scheme consisting of four equally spaced excitations with alternating x and y axis and α = 30 and γ = 175.
Both the " α x − α y " scheme (with large flip angle) and the " α x − γ y − α y − γ x " scheme preserve topology and also limit the effect of β.
Due to the alternating excitation axis, these composite propagators enforce a 3D trajectory of µ (see Eq. 11) for any β . In order to maintain the 3D trajectory for every possible β value, a minimum two RF excitations are necessary with different excitation phases. The large excitation angle for the " α x − α y " scheme may not be practical due to high demand on RF peak power when a sharp slice profile is used. A split of flip angles into two excitations can be realized in the " α x − γ y − α y − γ x " scheme: one of sharp slice profile, α , and one with no or very weak slice select gradients, γ.
Numerical optimization-choice of parameters in the selected scheme. The schemes were compared and the parameters of the schemes optimized with a Cramér-Rao lower bound (CRLB) analysis 15 .
Specifically, we evaluated the CRLB for all ten unique k = 4 schemes differing in the axis around which the RF pulses are played out, using two alternating flip angles α ∈ [20 • , 180 • ], γ ∈ [10 • , 190 • ] , with steps of 5 • . The simulation used the actual RF profiles as played out on the MR scanner scaled to the nominal flip angle.
For all evaluated settings, the coefficient of variation (CV) of the T 1 and T 2 was evaluated as the division of the square root of the CRLB by the nominal T 1 (800 and 1100 ms) or T 2 (50 and 200 ms) values for β ∈ [0, 2π] . These T 1 and T 2 values approximately span most of the commonly quantified tissues 24 .
It was observed that for the schemes with low CV, the CV depended only weakly on B 0 . However, for some schemes strong dependence on FA was observed close to the minimum CV. Avoiding those situations, overall the " α x − γ y − α y − γ x " scheme with α = 30 o and γ = 175 o was considered to provide the best compromise over the range of T 1 , T 2 , and B 0 .
Hence this scheme is the building block of the new sequence that we named "Multi-Phase balanced, non-Steady State Free Precession" (MP-b-nSSFP).
Acquisitions. The IRB (ethics committee name: "Medische Ethische Toetsings Commissie Erasmus MC", https:// www. erasm usmc. nl/ nl-nl/ pages/ metc) approved the in-vivo study (protocol 2014-096) and the acquisition was carried out after obtaining the informed consent from the volunteers. All experiments were performed in accordance with relevant guidelines and regulations. Only one representative volunteer is presented in this work to show the outcomes obtained. The acquisitions were performed on a 1.5 T clinical scanner (GE Optima MR450w, General Electric Medical Systems, Waukesha, WI) with a 16 channel head and neck coil.
The reference map were obtained with a multiparametric method that currently is a product for different vendors (MAGIC-Magnetic Resonance Image Compilation-for GE scanners). It is based on QRAPMASTER 25 that uses a multi-echo acquisition of a saturation-recovery sequence combined with a Fast Spin-Echo (FSE) as readout to obtain quantitative maps of T 1 , T 2 , and Proton Density (PD). Once the quantitative maps have been obtained, weighted images are synthesized from these maps. The acquisition was performed in axial orientation (AC-PC) with a TR of 4.7, FOV of 31 cm, and voxel-resolution of 1.2 × 1.2 × 5.0 mm. The total acquisition time of 20 slices covering all the brain was 5 min and 34 s. The phantom was scanned with the same protocol.
Additionally, reference B 1 and B 0 maps were obtained for the phantom, in order to compare them from those derived from the method poposed in this work. The B 1 map is a two-dimensional gradient echo based pulse sequence as described in 26 . B 0 map is based on 2D GRE pulse sequence repeated two times with different TE as described in 27,28 .
The images with our sequence-MP-b-nSSFP-were acquired with a Field of View 24 cm, reconstruction matrix 256 × 256, slice thickness 5 mm and 64 rewinded spiral out arms with 512 samples per arm for complete data collection. The images were reconstructed by density compensated non-uniform fast Fourier transform (FFT).
For the proposed MP-b-bSSFP, we acquired 25 repeats of the α x − γ y − α y − γ x scheme (Fig. 3e) α = 30º, γ = 175º, four readouts per block and T R = 120 ms per block with a delay of 3 s between blocks that acquire different spiral arms. The total acquisition time was 6 m and 40 s in this proof of concept study in which no acceleration was used. Due to radiofrequency specific absorption rate (SAR) and time restrictions, the γ pulses were not slice selective and instead surrounded by crushers to suppress out-of-slice signals.
The maps were obtained by for each voxel fitting where S is a vector with the complex valued signal of all 100 echoes, Bloch(θ) is a single species bloch simulation using hard RF pulses, E1 = e −TR/T 1 , E2 = e −TR/T 2 model the longitudinal, respectively transversal magnetization decay, β the offresonance induced phase evolution, B 1 a scaling factor for the RF pulses and S0 is proportional to proton density and contains the transmit-receive phase. After fitting, the T 1 and T 2 were computed from E1 and E2. The non-linear optimization was started from the best 5 out of 1000 candidate points, pseudo-randomly selected in the range E1 ∈ [0.91] , E2 ∈ [0.31] , β ∈ [−ππ] , B 1 ∈ [0.51.5] , with linear least squares solution for S0 . The non-linear optimization was performed with a custom trust region quasi Newton method implemented in MATLAB and the final fit with the lowest cost was returned 29 . The entire fitting procedure took 0.56s/voxel at our workstation (i7-8700 CPU).
An evaluation of the applicability of the model to the acquired data in the presence of intra-voxel B0 dispersion is presented in the supplementary material. Figure 4 shows the typical banding artifact present in balanced sequence and how the theoretical model (see Eq. 1) closely matches the data points despite the presence of external field ( B 0 ) inhomogeneities. www.nature.com/scientificreports/ Figure 5 shows images and different time points of the signal evolution and the signal evolution for three different tissues using as propagator the " α x − γ y − α y − γ x "scheme in a balanced pulse sequence. The banding artifact is gone and the signal evolution is also described by the same theoretical model (see Eq. 1).
Next figure shows the T1 and T2 maps from the MP-b-nSSFP and MAGIC (Fig. 6). Figure 7 shows a Bland-Altman comparison between MAGIC and MP-b-nSSFP.
In-vivo parametric and synthetic maps. Figure 8 shows the results from the MP-b-nSSFP in-vivo scan including the HO maps explained above. Additionally, Fig. 6 shows synthetic weighted images 25 . The T 1 -weighted image was simulated with TE = 20 ms and TR = 300 ms. The T 2 -weighted image was simulated with TE = 120 ms and TR = 4500 ms. The T 2 -FLAIR images were simulated with TE = 120 ms, TR = 15,000 ms and inversion time (TI) equal to 3000 ms. In order to provide a fair comparison we have used with the same protocol described above. Figure 9 shows the maps as acquired with MAGIC. Table 1 summarizes the estimated values of T 1 and T 2 in gray and white matter with MP-b-nSSFP and MAGIC and reference values from DESPOT1 and DESPOT2 30 , PLANET 31 and QRAPMASTER 25 respectively.

Discussion
In this work, we conceptually addressed several important aspects for enabling quantitative mapping from the fast transient response in a balanced pulse sequence. With this we extend previous work in this domain 8,14 to allow composite propagators. It deviates from MR fingerprinting approaches 10,11 in which experimental parameters are varied along the transient evolution and relaxation parameters are estimated by dictionary matching.  This comprehensively provides simultaneous mapping as well as maps of the experimental conditions. We proposed a train of acquisition blocks without variations along the train and without suppression of any part of the signal. We showed that a simple and very compact description can be provided for such a sequence. To our knowledge, the resulting real-valued expression has not been previously published.
Additionally, there were two separate aspects that we challenged. One was the impracticality or even impossibility to describe the signal in an analytical form in the transient state, especially when no part of the signal is suppressed (e.g. spoiling). The second aspect was the unfeasibility of using a single species approximation of complex and heterogeneous voxels for the transient evolution of a balanced sequence.
Moreover, we simultaneously addressed a long-standing problem of the balanced Steady State Free Precession technique: banding artefacts. By virtue of the conceptual design of the sequence (i.e.: 3D trajectory maintained by the composite propagator) every point in the image plane bears the same signal evolution characteristics. Along the transient evolution the phase dispersion is limited (see Fig. 3d,e) as the oscillation frequency ( ϕ , Eqs. 11, 12) is nearly the same for all off-resonance frequencies ( β) . Consequently no (propagating) banding waves in the image space will appear at any time point of the transient response (compare Fig. 4 vs Fig. 5), achieving the desired artefact-free images from balanced MR sequences. This could also be interpreted as a result of the γ pulse which acts as a refocusing pulse, however, the sequence is still sensitive to B 0 inhomogeneities across the image because the model contains and estimates the intra-voxel phase accrual. To our knowledge, this is a novel achievement that is different from the typical phase cycling in different runs 32,33 .
The resulting signal evolution demonstrates the validity of the single species description despite experimental imperfections.
It can be observed in the results from a standardized phantom (Fig. 7) that most of the results for T 2 are in good agreement with MAGIC. However, the results for T 1 are underestimated. Possibly, this is due to slice profile imperfections [34][35][36][37] and finite echo train length as well. www.nature.com/scientificreports/ In "In-vivo" experiments, T 2 results are also in good agreement wih MAGIC. For T 1 , there is a substantial underestimation as well. There are some additional intra-voxel physiological conditions that could be hindering the accuracy in T 1 : intrinsic intra-voxel assymetries [38][39][40][41] , water molecules exchanges 42 , diffusion 43 and magnetization transfer 31,44,45 . If validated and specific enough, such sensitivity to these other phenomena could be exploited for obtaining physiological information such as myelin content. The supplementary material demonstrates that in the current implementation intra-voxel B 0 dispersion and off-resonance can bias the T 1 and T 2 estimation. However, still the effect of B 0 dispersion is limited as the estimated T 2 value is substantially above the T * 2 value that would be obtained in, for example, a gradient echo experiment. Addressing these T 1 and T 2 biases is a topic of further research.
Our method relies on the topological interpretation of an algebraic object (i.e. the propagator: a transformation between the parameter space and the eigenvalue space). With the proposed propagators preserving the topology, the analytical single species description remains valid and allows simultaneous multi-parametric mapping.
It allows the separation of conceptual details of the propagator design and the numerical optimization of the nominal experimental parameters. Regarding the propagator, our method allows single species description even on the presence of B 0 inhomogeneities based on the range of the appropriate eigenvalues without excluding those complex (unlike Asslander's proposal 15 ). Furthermore, the 3D trajectory of μ is maintained at any value of B 0 in order to avoid loss of information (none of the normal mode vectors vanish). Moreover, the mapping between parameter space and eigenvalue space preserves topology in order to distinguish between species by their signal evolution. To demonstrate the feasibility on actual scanners, we have selected one of the propagators that fulfills these requirements. In addition to the analytical approach, a numerical optimization of parameters such as flip angles and phases is carried out for the best estimation of T1, T2, PD, B1, B0. So, our novel approach not only enables the estimation of intrinsic parameters, but from the same data and estimation process, one can derive the imperfections of the experimental setup. The macroscopic confounding factors (experimental imperfections such as B 0 inhomogeneities and B 1 inaccuracies) are not simply suppressed or excluded in our acquisition, rather they are included in the theoretical description as parameters to be determined (see results "d" and "e" in Fig. 6).  www.nature.com/scientificreports/ To our knowledge, no other parametric imaging methods provide estimation of this set of intrinsic and experimental parameters simultaneously based on an analytical solution describing the transient response of the magnetization with using dictionaries.The concurrent estimation of intrinsic and experimental imperfections based on the proposed comprehensive analytical model makes this new technique less affected by common system imperfections and could allow for the development of less demanding MR scanners. Based on the robustness of the synchronized estimation of experimental imperfections and parametric maps, synthetic MR images for different clinically relevant contrasts can also be reconstructed without relevant artifacts (see results "g", "h" and "i" in Fig. 8). As a consequence of these features, our method could contribute to the standardization of MR imaging based on simultaneous multi-parametric mapping and synthetic weighted MR.

Conclusions
Building on signal-based MR, we provide a complete and comprehensive analytical expression for the signal evolution of a balanced sequence. We extend this solution from the continuous Bloch equation to a discrete description of an entire imaging sequence. The analytical expression relies on a simple, single species model of an imaging voxel which is shown to be appropriate despite the heterogeneity of voxels in-vivo. We demonstrate the importance of an analytical approach in the design of the sequence propagator. This theoretical model could be fitted to experimental data without requiring a dictionary. We simultaneously derive parametric maps of the intrinsic properties (T 1 , T 2 , PD) as well as from imperfections of the experimental parameters (B 0 , B 1 ). We demonstrate the feasibility of our method on clinical MRI scanners for in-vivo brain scans.